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Abstract 

The concept of robustness of regulatory networks has been closely related to the nature of the in- 
teractions among genes, and the capability of pattern maintenance or reproducibility. Defining this ro- 
bustness property is a challenging task, but mathematical models have often associated it to the volume 
of the space of admissible parameters. Not only the volume of the space but also its topology and ge- 
ometry contain information on essential aspects of the network, including feasible pathways, switching 
between two parallel pathways or distinct/disconnected active regions of parameters. A general method 
is presented here to characterize the space of admissible parameters, by writing it as a semi-algebraic set, 
and then theoretically analyzing its topology and geometry, as well as volume. This method provides a 
more objective and complete measure of the robustness of a developmental module. As an illustration, 
the segment polarity gene network is analyzed. 

1 Introduction 

For biological networks, the concept of robustness often expresses the idea that the system's regulatory 
functions should operate correctly under a variety of situations. The network should respond appropriately 
to various stimulii and recognize meaningful ones (either harmful or favorable), but it should also ignore 
small (not meaningful) variations in the environment as well as inescapable fluctuations in the abundances of 
biomolecules involved in the network ||Tl|2l|3J. One might even speculate that if the networks malfunctions 
easily as a result of mutations then it has low chance of being selected by evolution. In that case one might 
expect a certain degree of mutational robustness fT,"?]. 

While it is difficult to define this robustness property in a precise form, it has been associated to the space 
of admissible kinetic parameters, its volume [3], and the effect of paramater perturbations on the qualitative 
behavior of the system ||lJ|2l. Some methods for parameter sensitivity have been developed 15101, based 
essentially on derivatives of variables or fluxes with respect to the system's parameters. The volume of 
the parameter space can be used as an indication of "how many" parameter combinations are possible, and 
these are related to the ability of the network to work under a variety of situations. For instance, parameters 
may range through different orders of magnitude, representing very different environments. However, size 
is often not a reliable measure for robustness; other quantities, such as shape, play a much more important 
role, as illustrated in Fig. [T] Analysis of the shape or geometry of the admissible parameter set gives an 
indication not only of its size, but also how far perturbations around each parameter disrupt the network. 
A robust biological network will admit small fluctuations in its parameters without changing its qualitative 
behavior. So, a robust network will be associated to a system whose parameter set has few "narrow pieces" 
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Figure 1: The role of geometry and topology in robustness. Regions (a) and (b) have the same volume, 
but (b) is less robust: the same perturbation leads out of the space. Regions (c) and (d) also have the same 
volume, but (d) is not a simply connected set, hence less robust. 

and "sharp corners". In such sets, reasonable parameter fluctuations may occur without leaving the set, 
hence maintaining the network's qualitative behavior (compare Fig. [T] (a) and (b)). We can formalize a 
measure of robustness that is related to having low rate of exit from the region under random walk H . The 
rate of first exit is intrinsically connected to the geometry of the region and is particularly sensitive to narrow 
directions and not just the overall volume. 

To illustrate the importance of parameter space geometry, and the insight it brings to understanding the 
network, the model of the segment polarity network developed by von Dassow and collaborators 131 will be 
analyzed. The segment polarity network is part of a cascade of gene families responsible for generating the 
segmentation of the fruit fly embryo fT|. Genes in earlier stages are transiently expressed, but the segment 
polarity genes maintain a stable pattern for about three hours. It has been suggested that the segment polarity 
genes constitute a robust developmental module, capable of autonomously reproducing the same behavior or 
generating the same gene expression pattern, in response to transient inputs E [H 13 . This robustness would 
be due to the nature of interactions among genes, rather than the kinetic parameters of the reactions. The 
model [3 1 describes the interactions among the principal segment polarity genes, is continuous, and involves 
cell-to-cell communications and around 50 parameters which are essentially unknown. The authors of S 
explored the model by randomly choosing 240,000 parameter sets out of which about 1,192 (or 0.5%) sets 
were consistent with the generation (at steady state) of the wild type pattern. To explore the robustness of the 
network as a property of its interactions, Albert and Othmer |9| developed a Boolean model of the segment 
polarity network, a discrete logical model where each species has only two states (0 or 1; "OFF" or "ON"), 
but no kinetic parameters need to be defined. This Boolean model is amenable to various methods for 
systematic robustness analysis lfT0l[TTl[T2l . Ingolia [8 | focused on the properties of the (slightly changed) 
model [3 | in individual cells, such as bistability, and extrapolated necessary conditions on parameters to the 
full intercellular model. 

We propose a different approach, that retains the information contained on the kinetic parameters but 
reduces the model to a logical form with various possible ON levels and species-dependent activation param- 
eters. The admissible set of parameters of the model [3] is analyzed by constructing a cylindrical algebraic 
decomposition. Among other conclusions, our analysis completely explains the two "missing links" in von 
Dassow et. al. original model, namely: why the segment polarity pattern can not be recovered without the 
negative regulation of engrailed by Cubitus repressor protein, and why the autocatalytic wingless activation 
pathway vastly increases the network robustness. The present approach shows that, in contrast to volume 
only estimates, the topology and geometry of this set provide reliable quantitative measures of robustness of 
a system. 
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2 Steady states define the feasible parameter space 



Previous studies 111 [H have tested the parameter space by randomly choosing sets of parameters and simu- 
lating the continuous model. If the corresponding trajectory reaches a steady state, and if this steady state is 
compatible with the experimentally observed wild type gene pattern, then the given set of parameters is said 
to be a "solution" to the modeling problem. 

A more efficient and complete study of the parameter space can be devised, by first solving the algebraic 
equations of the model at steady state, and writing the steady state solutions as a function of the parameters. 
On the other hand, the steady state solutions are known - the set of elements representing the wild type 
pattern is denoted by W - so, one can then look for parameters that yield this pattern. Since many sets of 
parameters may be expected to yield the wild type pattern, this procedure provides a family of conditions 
defining regions of "good"or feasible parameters "p" for wild-type steady states x G W. 

The von Dassow et. al. model Before proceeding, recall that the model (Appendix |B]| describes the 
concentrations of various mRNAs and proteins in a four cell parasegment of the fly embryo, subject to 
periodic boundary conditions (see also Fig. |2]). Here, each cell is assumed to have a square shape, with 
four faces (see Appendix [E]|. We next very briefly recall the species involved. There are nine species with 
homogeneous concentration throughout each cell: engrailed mRNA and protein {en and EN), wingless 
mRNA and (internal) protein (wg and IWG), patched mRNA (ptc), cubitus mRNA, active and repressor 
proteins {ci, CI, and CN), and hedgehog mRNA (hh). Each of these species has a distinct concentration in 
each cell (Xi, i = 1, ... ,4). In addition, there are three other species whose concentration varies in each 
of the four cell faces: external wingless protein (EWG), patched protein (PTC) and hedgehog protein (HH). 
For each of these species, the concentration in cell i at face j is denoted Xij, i = I, ... ,4, j = 1, . . . , 4. 
Thus, overall there are: n = 9x4 + 3x4x4 = 84 variables. Throughout the paper, the following notation 
will be used (prime denotes transpose): 

X = {Xi,X2,X3,X4y, forX E {en,EN,wg,lWG,ptc,ci,G,CN,hh}. 

and 

X = Xi,2, Xi,3, Xi,4, X2,i, ^4,4)', for X e {EWG, PTC, HH}. 

The total vector of concentrations is: 

X = {en, EN', wg' , IWG', EWG', ptc, PTC', ci' , Cl' , CN', hh' , HH'). 

Set of feasible parameters In general, the problem can be formulated mathematically by writing a set 
of equations dependent on the vector of species concentrations (x G K>o) the parameter vector (p € 
M>q), together with a set of outputs (y € ffi>o, the available gene expression levels). Introduce functions 
f~.Wi,Q X W and h : M^q y CWI^q, where M>o = {x eR: Xi>0, for afl i}, and consider 

the system with outputs 

f = fi^^P) (1) 
y = h{x) (2) 

where the function h{x) could be, for instance, a vector listing the concentration of wingless, engrailed, 
hedgehog and cubitus, four of the segment polarity mRNAs which have been experimentally measured. Or, 
in other words, y is "the phenotype corresponding to the genotype x". The wild-type gene expression output 
set can be defined as: 

y-^^ = {yey: y = h{x), X G W}. 
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The problem of characterizing the sets of feasible parameters is then reduced to finding all possible param- 
eter vectors p which lead the system to have an output in 3^*^, at steady state. This will be the set of "good" 
parameters: 

G = { p G M>o : 3x s.t. f{x,p) = and h{x) G 3^*^ }. (3) 



Large Hill coefficients A straightforward approach would be to solve the original system at steady state, 
obtain expressions for a; G W in terms of p, and compare these expressions to the outputs in J^*^: 

f{x,p) =0 4^ x = F{p) and F{p) G W <^ p G G. 

A possible drawback of this method is that explicit solutions x = F{p) for the original system and then 
explicit formulas for G may not be easy to compute. On the other hand, many of the equations in the 
model [3J involve terms of the form (see also Appendix [B|): 

(p{X,K,iy] 



+ X"' 



meaning that the function (f) is active (ON), if species X is above a certain threshold n. The exponent u, 
also known as the Hill coefficient, characterizes the steepness of an OFF/ON transition. For large enough 
exponents, this saturation function becomes very steep, and cj) becomes practically insensitive to the actual 
value of u. As found in [13 |. coefficients v must indeed be quite large for the network to achieve robustness: 
namely in the interval [5.0, 10.0]. This is also the basis of the typical on/off logical interpretation of gene 
expression. Any such term 4'{X, k, u), for large u, may thus be replaced by a step function with two levels 
(0 or 1): 



e{x 



0, X < K 

1, X > K. 



Thus, when v is large: 

lim (j){X,K,iy) = e{X - k), lim ip{X,K,i^) = l-e{X - k) =9{k- X). (4) 
A composite function of (p and ip also frequently appears in the continuous equations (Appendix [B]): 

This function can be simplified in terms of step functions to: 

e{Xae{Kb - Xb) - Ka) = e{Xa - Ka)e{Kb - Xfe) 

since 



Xfe > Kfe ^ e{Kb - Xfe) = ^ 0{Xae{Kb - X^) - Ka) = 9{-Ka) = 0, 
Xh<Kb^ e{Kb -Xb) = l^ e{Xa9{Kb - Xb) - ^a) = 0{Xa - Ka) . 

As an example, consider the equation governing engrailed from the original model which can be found 
in |l3l [131 (or in Appendix |B]|. In this model the concentration of engrailed in cell i {eni), is positively 
regulated by external Wingless protein (EWGj) and negatively regulated by Cubitus repressor protein (CNj) 
concentrations (further notation is found in Appendix [A]): 

dcfi ' 1 

-T^ = TP ((/)(EWGi'0(CNi, Kcnen, t'cN.J, ^wG™ , i^wG™) " eni) . 
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For large exponents v, this simplifies to the equation: 
derij 1 



dt H„, 



(0(EWGi - Kwg.,)^(kcn™ - CNi) - erii) 



To analyticaly study the space of feasible parameters for the segment polarity network model f3l , we will 
thus consider that all exponents v are large, and apply method (|4]) to simplify the original system of equa- 
tions. The von Dassow et. al. model is then characterized by equations ([3T])-([42]) (Appendix [C]). The pa- 
rameters are as in Q, except Ti and Ui, which represent the maximal values of ptc and ci (respectively), 
in each cell. These take values in the interval [0, 1] and generalize the possible ON values of ptc and ci (to 
be discussed later). In addition, as discussed, the system is assumed to be at steady state, in which case the 
gene expression pattern must satisfy: 

erii = 9{EWGi - Kwg™)^(«^cn<.„ - CNj). 



Applying Q and then solving the system at steady state yields the set of algebraic equations (|43|)-(|54|), which 
characterize the gene expression pattern of the segment polarity network according to the von Dassow et. 
al. model. 

Maximal (ON) expression levels While some of the species have a normalized maximal expression level 
(to 1), such as en or hh, other species may be more generally allowed to have any positive value (namely, 
ptc and ci). These maximal expression levels are also treated as parameters. When using (HI) to simplify 



the patched equation (24i to (36 1, we have generalized the equation and added distinct maximal levels 
of expression in each cell, given by Ti {i = 1, ... ,4). This allows a more accurate representation of 
experimental data, which shows that patched is strongly expressed in every second and fourth cells, weakly 
expressed in every first cell, and not expressed in every third cell (see |i3J for more discussion). Thus we will 
consider Ti<T2=T4: 

P<2,3,4 = (^l'^2,0,T2)'. (5) 

A similar generalization was made to deal with the activation of cubitus interruptus. In von Dassow et. 
al. model, this is due to some external parameters Bi (not governed by a dynamical equation), with a 
corresponding activity threshold Kb^,. However, for more generality, and to allow distinct maximal levels 



of expression in each cell, we have replaced each of the terms 9{Bi — Kbc,) in (38 1 by a parameter Ui, 



i = 1, . . . ,4 (50l. Furthermore, in characterizing the set of feasible parameters, it will become clear that 
allowing distinct Ui enlarges the space of possible parameters, by introducing the four regions Gc,i to Gc.iv. 
Thus the steady state values for the cubitus mRNA are: 

<,2,3,4 = (f^l,f^2,0,C/4)'. (6) 

Asymmetry in cubitus expression (i.e., distinct values Ui) could be due, for instance, to some of the pair 
rule genes. Sloppy paired, or a combination of Runt and Factor X, regulate the transition from pair rule to 
segment polarity genes expression, and induce asymmetric anterior/posterior parasegment expression [14|. 
Finally, note that the maximal expression levels of wg are expressed in terms of the parameters acug and 



OwGug- From equation (45 I, there are three possible combinations of the step functions, each leading to a 



different value for wg2- These three possibilities are: 

''^Cl — T~~ ) 'ii'cijWG — T~~ ~~ ) W^Q — ——- , 

i + OLciKg i + Olcwg ~r ttwGii'g ^ ~r OwGivg 

and each reflects a different pathway for wingless activation. Indeed, wingless can be activated by Cubitus 
only (in which case the maximal amplitude is given by Wa), by both Cubitus and Wingless (u'ci.wg), or by 
Wingless only {w^a). 
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Outputs The next question concerns the choice of an appropriate output function. The gene expression 
patterns for engrailed, wingless, hedgehog, cubitus, and patched are among the most well documented, so 
we will consider the output function h : M>q 



!>20 . 



y = Kx) 



(KXx)\ 

ha{x) 



(V) 



At steady state, both en and hh are expressed in every third cell lITSll . which translates into 

h,„{x) = (0, 0, 1,0)', hi,i,{x) = (0, 0, 1,0)', for x G W. 



(8) 



Further experimental observations show that cubitus is expressed in all but the third cell [16], and patched 
is strongly expressed in every second and fourth |[T5l . but more weakly expressed in every first cell. So: 



K,{x) = {Ui, U2, 0, C/4)', K,{x) = (Ti, Ta, 0, Ta)', for xeW. 
Finally, wingless (wg) is only expressed in every second cell lITSll . to the left of en, that is: 

h„^{x) = (0, w, 0, 0)', for xeW. 
To summarize, in this example, the set of output values at steady state is: 

3^*^ = {((o,o,i,o),(o,ti;,o,o),(ri,r2,o,r2),([/i,c/2,o,f/4),(o,o,i,o))' : 

w, Ti, T2, Ui, U2, f/4 > 0, Ti < T2} . 



(9) 



(10) 



(11) 



The first result to be noted is that there is a unique steady state x = x{p) G W for each set of parameters p: 



Theorem 1. Let f be the function M>q x ]R>q M" given by {31 \-{42 \, and h be the function M>q y 
given by Define G as in (jjj). Then, there exists a function F : G M>g such that, for each p £ M>q 



and each x £ 



f{x,p) = and h{x) £ imply x = F{p). 



Proof. Pick any p £ G, and an x G M" q satisfying f{x,p) = and h{x) £ 3^*^. The equations f{x, p) = 
can be simplified to yield (43 l-(54i. We must check that these equations are all consistent and admit only 
one solution. Since all mRNAs en, wg, ptc, ci, and hh are provided by h{x), we must solve for the proteins, 
and then substitute these back into the equations for the mRNAs, to check consistency. 

The Engrailed protein is straigthforward: EN = en. We start by solving for EWG, with wg = 
(0, w, 0, 0)' as given. First note that the matrix M is diagonally dominant, by adding up the entries in 
any column: 



(^iwG + '^endo + + 2r2,Af) + '^tlm + rM + = -H.J. 



1 



1 + i^n 



(12) 



which is always a negative quantity. By Gersgorin's Theorem, all eigenvalues of M are contained in the 
disk centered at —d + h with radius 2rLM + i^m + 3/i, therefore all eigenvalues have negative real parts. 
Thus, the matrix M is symmetric and negative definite, and since the right-hand-side vector in ( 47 ) is also 
non-positive, there is a unique solution 



EWG 



1 T 

' exo 

4 1+ Hiv/arc: 



M"^ wg 
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which is real and positive, /or each set of parameters p. Once we have EWG, we can immediately solve for 



IWGfrom(46l. 



The solution for PTC and HH can also be exactly and uniquely computed from ( |49] ) and ( [54] ), for any 
output ptc = (Ti, T2, 0, T2)' (this calculation is shown Appendix|F]l. 



Finally, one can now straightforwardly and uniquely compute the values of CI and CN, from ( 5 1 1 
and ^52), and the values of ci and PTC. 



The last step is the substitution of EN, EWG, IWG, PTC, HH, CI and CN back into the equations for 



the mRNAs (43 1, (45 1, (48 1, (50 1, and (53 1. But, since p G G, by definition we are guaranteed that these 



equalities are indeed satisfied. I 

Missing link: engrailed regulation by Cubitus repressor A second result from our model formulation 
is the explanation of a "missing link" in a first version of the model proposed by von Dassow et. al. Q. 
In this first version, engrailed was regulated only by EWG, and no feasible pai^ameter sets were found. 
Indeed, below (Theorem [2]l we prove that, /or any set of parameters, the mechanism for wingless regulation 
generates a strong symmetry in the steady state expression of external Wingless. This symmetry effectively 
prevents any asymmetry arising in en due to EWG only. 

Theorem 2. Let w > and assume wg'^^ = (0, w, 0, 0)'. Then, at steady state: 

EWGl^ < EWG^^ = EWC^^ < EWG^. (13) 



The proof is based on a sequence of algebraic calculations, and is shown in AppendixjE] Now, consider 
the steady state equation for engrailed, when no dependence on CN is assumed: 

enr = d{EWGr - a^wg.,) 

Compare to the output ([8]): 

K„{x) = (0,0,1,0)'. 

Then, from the definition of 0, for consistency in our model it is necessary that: 

EWG7^<Kwg™, fori = 1,2,4 
EWG^ > KwG.„. 



However, by (13 1, the inequalities for i = 1,2 and i = 3 are incompatible. This means that, due to the 
symmetry in Wingless distribution, such a simple regulation of en can never lead to the segment polarity 
pattern. Thus engrailed requires regulation by some other factor, in this case repression by the Cubitus 



protein (CN), as in (43 1. In order to obtain repression of en in the first and second cells, one can now ask: 



cNr > 

"'CNc/i 

EWGY > Kwg™ and CN^^ < Kcn™ 
EWCr < KwG„, orCN7 > KcN., 

that is, CN is responsible for repression in both the first and second cells. This means that, at steady state, 
CN must be expressed in both the first and second cells. This in turn requires the presence of Patched protein 
in both the first and second cells. On the other hand, from Appendix [Fj we know that a steady state x & W 
with h{x) G implies ptc^^ = PTC^, and also PTC^^ = PTC^. This can be stated as: 

Lemma 2.1. Consider system ([T]) and assume that, at steady state, the output set is J^^^. Then ptc^^ = 
PTCr / 0. If pte*'^ = (Ti, Ta, 0, Ta)' with Ti < T2, then PTC^ = Ti and PTC^^ = PTC^ > 0. I 
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While patched expression is typically weaker in the first than in second and fourth cells (see [3]), this 
shows that it is nevertheless necessary, that is, the segment polarity gene pattern obtains only when T\ > 0. 
The discussion on CN leads to the following conclusion: 



Lemma 2.2. Consider system ([T]) and assume that, at steady state, the output set is y 

(Ti, Ts, 0, Tz)' with Ti < Tg. Then VTCl\ > Kptcci and 



Let pfc" 



cir 



1 



1 + HriCr 



CNf 



HriCr 



1 + HriCr 



1,2,4, 



(14) 



and Cir 



CNT 



0. 



Proof. Theorem|2]and the subsequent discussion shows that CN72 / is needed. From (52 1, this can only 



be achieved by asking PTC{^2 > ^ptcci- By Lemma 2.1 it also holds that PTC2 



means that both ( |51| ) and (52 1 can be simplified to ( 14 1, at steady state. On the third cell, CI 
because the output is zero. 



PTCJ^ 



^ '^PTcci- This 




CN^" 



3 A cylindrical algebraic decomposition of the parameter space 

The algebraic equations f{x,p) = together with h{x) G 3^*^ are a representation of the set of good 
parameters G, though not providing as yet explicit conditions on p. An explicit characterization of the 
parameters p may be obtained by calculating a cylindrical algebraic decomposition (CAD) of G: this is a 
special type of representation of G as a finite union of disjoint connected components. A CAD will provide a 
hierarchy of inequalities on pi, p2,. .., Pr, from which the volume of G, as well as its geometry and topology, 
may be deduced. 

Computing the cylindrical algebraic decomposition of a semi-algebraic set is a complex problem, but 
various standard algorithms are available ifTTlfTSl . Several software packages have been developed, for in- 
stance QEPAD [19], (based in [20]) and in Mathematica [21 J. See also [22J for an overview of available 
software, current applications, and many other related references. Common applications of CADs include 
computation of the controllable or reachabable sets in hybrid systems ||23ll . Constructing a CAD involves 
the use of symbolic computation and, while various improvements have been achieved, it still is a time con- 
suming problem. For instance, the estimated maximum time for the algorithm [ 17J is dominated by "2^ ", 
where N is the length of the input formula and < A; < 8. Fortunately, in view of these computational 
complexity difficulties, in the present example it is relatively easy to directly compute a CAD without using 
general methods, and we will do so. 

For equations (43 1-([54)), subject to ([8])-( 10 1, a cylindrical algebraic decomposition can be constructed 



in which several parameters (Table |6]l are free to take any values (within physiological restrictions only). 
At the next level, parameters in Table [7] have constraints which depend only on those parameters given in 
Table [6] The last level is formed by the parameters in Table [8] whose constraints depend on parameters from 
both previous levels (Tables |6] and |7]), thus defining a polyhedron. 

Following the model of von Dassow et. al., there are two possible parallel pathways for wingless acti- 
vation: either by the Cubitus interruptus protein (CI), or through auto-activation; both pathways could be 
simultaneously activating wingless production. Since the activation constants aci„g and a^^,G„■g, are free pa- 
rameters, in each of the three cases wg2^ will have a different ON level (respectively, Wa, Wwg> or Wci,wg)- 
Computation of EWG and IWG depends on wg*^, so each of these three cases must be separately analyzed 
for feasibility. For both pathways, exact analytic computation of PTCjj and HHjj = 1, . . . , 4) is also 
carried out (see Appendix [F]). Several disconnected regions of parameters will be defined by the levels of 
cubitus, C/i 2 4- 
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Five disconnected regions When only CI and CN regulate wingless expression, it is easy to see from ( 45 1, ( 10 1 
and (fT4]l that: 



1 + HciCci "^'"^ * 1 + HciCc 



Uil—FT-TT- < i^ci.s or Ui > Kcnh^; and IWGi < Kwg„s (15) 



for z = 1, 3, 4, and 



1 + HciCci 1 + HciCci 



U^l—rr-TT- > /^ciHi; and < Kcn..^ and IWG2 < Kwo^s • (16) 



From observation of ( [T5| ), ( [T6| ) it is clear that the situations U2 = Ui or U2 = U4 are not well defined, since 
contradictory constraints are imposed on Kci„.g and KcNug- So, the regions of parameters satisfying U2 = Ui 
or U2 = U4 are not feasible. This divides the set G into at least four disconnected components, divided by 
the hyperplanes U2 = Ui or U2 = f/4 (Gci, Gen, Gc.m, and Gc.iv in Fig-[3]l. A similar argument holds for 
the case when both pathways contribute to activation of wingless on the second cell. The four disconnected 
regions of parameters are identified in Table [8] 

Finally, the third case (auto-activation pathway only), introduces a fifth component of G (Gauio). which 
must be disconnected from either of the previous four components. This is clear, by contrasting the necessary 
conditions in the second cell for either case (compare to ([16])): 



^ ^ or [/ -^"^c 

1 + HciCci 1 + HciCci 



U2 , , ,r ^ < i^cin-s or U2——fr-pT- > '^cN„g and IWG2 > Kwong- (17) 



The five disconnected components are thus first defined by Ui, U2, and C/4, and then by Kq,,,, and KcNng- The 
projection on the (Kci„j,KcN„g,KwG.,g)-dimensions compares two of these components (Gen and Gauio), both 
polyhedrons (Fig.|4]l. 

The cylindrical algebraic decomposition is shown in detail in Appendix [Gj and summarized in Ta- 
bles 6|7|8 Each of the five components, G^ G {Gc.i, Gc u, Gcmj Gcjv; Gauio} is thus described by a hierarchy 
of sets of the form: 

Si^-j = {cL^i by) C M 

Si,-y = {{x,Xi) G M* : X e Si-i^y, ai^y{x) < Xi < Pi^-yix)} C (18) 

for i = 2, . . . , A^, where ai,^, f3i,y : Si^i^^ M>o and S]\j^j = Gj. It can be shown that each Gj is in fact 
topologically equivalent to the unitary open hypercube, and hence topologically trivial. 



Theorem 3. For each Gj, the set Sn = Sj\[^j, as obtained from {18), is homeomorphic to (0, 1)^. 

Proof. Pick any G^, and drop the subscript 7, for simplicity of notation. To argue by induction, note that 
the set Si is clearly homeomorphic to (0, 1). For i > 2, assume that Si-i is homeomorphic to (0, 1)*~^. 
Next, define the following continuous function: 

(pi : Si-i X (0, 1) Si-i X M, (pi{x, t) = {x, fi{x) + t {gi{x) - fi{x))). 

For each fixed x, ai{x) < ai{x) + t {(3i{x) — ai{x)) < Pi{x) for all t e (0, 1). Therefore, (fi maps into Si. 
On the other hand, ipi has an inverse function defined on Si and continuous, given by: 

: 5i ^ S'i-i X (0,1), Lf- [x,y)=[x, 



/5i(x) - ai{x] 

So Si is homeomorphic to Si-i x (0, 1), and therefore, by inductive hypothesis, to (0, 1)*, as we wanted to 
show. I 
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Relative volume and the second missing link Once the parameter set G is characterized by writing 
intervals for the various parameters in the form ( 18 1, it is very easy to compute the (relative) volumes of the 
disconnected components. Note that in each component only the intervals for Ui, U2, Ua, and Kcn<.„, Kwg™, 
«^ci,vg> KcN„g, KwGv,'g vary; constraints on the remaining parameters are common to all components. Following 
a Monte Carlo approach, the parameters in Tables 6]7 are chosen first, and then Kcwg, Kcn,,?, KwG«g from 
the unitary cube (all parameters are randomly chosen from uniform distributions in the given intervals). 
It is next checked whether the parameter set falls in any of the components Gc.i, Gcu, Gc,iu> Gcj\, Gp,^,^, 
or outside G. This method provides an estimate of the volumes of each disconnected component, when 
projected into the (KcNf.„,KwG«,>'«ciug>'icN.vgAwG«g) dimensions, as the fraction of parameter sets that fall into 
each component. The volume of this 5-dimensional cube occupied by feasible parameter sets is only about 
0.7%. As is illustrated by the polyhedrons in Fig. [4j component Gauio is much larger than the others - 
approximately 40 to 270 times larger. 



Table 1: Relative volumes of the five disconnected components. In component Gauio. only auto- activation 
leads to wingless expression. In components Gc,i, Gc.u, Gc.m, and Gcw, CI always activates wingless ex- 
pression. Total number of parameter sets generated: 1 x 10^. Number of feasible parameter sets: 70026. 



Component Volume 



Gc,i 


1.6 X 10" 


-4 


Gc,i[ 


0.25 X 10 


-4 


Gc,m 


0.86 X 10 


-4 


Gc,lV 


0.46 X 10 


-4 


c 

^Auto 


67 X 10- 


4 



The large difference observed between Gc.i-Gc.iv and Gai„o explains the second "missing link" in the first 
version of von Dassow et. al. model, namely the wingless autocatalytic activation. Note that the presence of 
this link greatly increases the total volume of the feasible parameter space: in fact the region Gauio is 95% of 
the total volume. 

Parameter tendencies As described above, the parameter space for the segment polarity network can be 
described by a CAD, a hierarchy of inequalities on the parameters where an interval is explicitly given for 
each parameter. At the base of this hierarchy, there is a first group of parameters whose intervals correspond 
simply to physiological values, as in Table[6] The intervals for the remaining parameters have bounds which 
depend on the parameters in the first group (Tables |7] and [8]l. In any case, one may ask how the parameters 
are distributed in their intervals, whether each parameter pi is more likely to attain high or low values more 
frequently, or whether a "tendency" for each parameter pi be identified. An answer to this question is 
obtained by randomly generating parameters in the full parameter space G, and computing the distribution 
of each parameter. Taking all the parameter sets generated to compute the relative volumes of the five 
disconnected components of G, and computing a histogram for each parameter, the result shown on Fig. [5] 
is obtained. As expected, many parameters have a uniform distribution, as their values do not influence 
the final outcome of the network in any particular way (for instance, most half-lives). Other parameters 
exhibit a marked tendency for higher (e.g., k,ch,„c), medium (e.g., KwGv,g) or lower (e.g., Kd,,,,) values. All the 
parameters that exhibit a marked tendency are listed in Table [2j and classified according to their function in 
the network: for instance, Kcn,,,^ represents the repression of ptc by CN, and therefore, high values of Kc^p,, 
correspond to a weak repression. 

A very similar analysis was performed by von Dassow and Odell |[T3]| . who also plotted the distribution 
of their family of feasible parameters to determine possible constraints for each parameter. Overall, our 
results agree very well with those of von Dassow and Odell: most tendencies found by these authors (see 
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Fig. 6 and Table 1 of (13]) are confirmed by our parameter space analysis. There are only five exceptions, 
where our analysis showed no tendency (compare columns 3 and 4 of Table |2]), suggesting that these five 
parameters can, in fact, take values in a larger set, implying that the parameter space is larger than estimated 
in lfT3l . From these exceptions, Kenci, Kew,i„ ^-cmi,, and re„joWG all belong to the group of parameters which 
can be freely chosen. The other parameter is Kd,,,, which depends on the disconnected regions, and again 
our analysis shows that this pair has no preferred tendency. 

A more detailed examination of the conditions on Kcurg and Kcn„s turns out to be very illuminating. 
First, note that Kq,,,, and KcNv.g define the five components, in the sense that distinct intervals for these two 
parameters are given in each component. Thus, it may be expected that the distribution of these parameters 
varies in each region (Table [3]). Indeed, by plotting the histograms for Kci„.j and Kcn„j for each region 
alone, we note that these show a marked tendency in components Gc.i — Gc.iv, for low Kci,,,^ and high KcN„g- 
In contrast, the distributions of Kd,,, and KcN«g for region Gauio alone show an opposite tendency. This is 
consistent with the fact that the volume of Gamo is about 95% of G and, therefore, it dominates the overall 
tendency. Note also that, in the four components Gc,i — Gc,w, it always holds that Kcihs < «^cNv.g, clearly 
in agreement with the tendency observed for our parameter sets. In component Gauio the parameters Kcih^ 
and KcNug must satisfy constraints that contradict those of Gc.i — Gc,iw, but not necessarily exactly opposite 
constraints (see Table [8]l. Thus more freedom results for the choice of Kciwg and KcNv,g in Gauio- The tendency 
of Kciwg and KcNug in Gci — Gc.iv is, however, the opposite of that observed by von Dassow and Odell, a 
fact that can be explained once again by the "second missing link". Indeed, since all feasible parameter 
sets in ||3] [131 were found only after adding the autocatalytic wingless activation link, it can be inferred that 
those parameters belong to region Gauio- We conclude that the parameter space is larger than estimated by 
von Dassow and Odell. 

Table 2: Comparison between the constraints identified by von Dassow and Odell |[T3l . and the exact con- 
straints given by the five regions defined above. Total number of parameters generated in G: 70026. 



Parameter 


Description 


Tendency 


Tendency 






(Ha, Table 1) 


(within G) 




WG activation of en 


Moderate 


Moderate 




CN repression of en 


Strong 


Strong 


I^WGwg 


WG autoactivation 


Moderate 


Moderate 


Kchvg 


CI activation of wg 


Weak 




I^CNwg 


CN repression of wg 


Strong 


Strong 


Kciptc 


CI activation of ptc 


Strong 


Strong 


I^CNplc 


CN repression of ptc 


Weak 


Weak 




EN repression of ci 


Moderate 




'^PTCCI 


PTC stimulation of CI cleavage 


Strong 


Strong 




EN activation of hh 


Weak 






CN repression of hh 


Strong 




Ccr 


Maximal cleavage rate of CI 


Rapid 


Rapid 




Half-life of intracellular WG 


Short 


Short 


'"endoWG 


Rate of WG endocytosis 


Slow 




'"exoWG 


Rate of WG exocytosis 


Moderately slow 


Moderately fast 


'^MxfeiWG 


Rate of WG cell-to-cell exchange 


Slow 


Slow 


O-WGtvg 


Maximal WG autocatalytic rate 




Moderately rapid 
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Table 3: Influence of the autocatalytic WG activation link in the parameter distribution. 



Parameter 


Tendency 


Tendency 


Tendency 




(113, Table 1) 


^Auto 


Gc.l, Gc,lly Gc_in Gciv 






(WG ^ wg) 


(WG ^ wg) 


KclHg 


Weak 


Weak 


Strong 




Strong 


(Moderately) Strong 


Weak 



4 Geometry and robustness 

The volume estimates for the parameter space regions give an idea of "how many" parameter combinations 
are possible. But volume alone is often not a reliable measure for robustness, as illustrated in Fig. [T] The 
shape or geometry of the parameter space regions also shows how far perturbations around each parameter 
will disrupt the network. Thus, parameter regions exhibiting "narrow" pieces or "sharp" comers indicate 
a lower level of robustness in the network. One way to explore the shape of a given multi-dimensional set 
is to consider a random point (p^) and follow a random walk in space (p^ = p^~^ + dp'', k = 1,2,...,), 
where each step has the same absolute length (\dp^\ = oq), but a random direction. Then record the number 
of steps needed for the point to exit the given set. Repeating this procedure for many points in the set, the 
probability that a point leaves the set after t steps can be computed. 

The random walk could be interpreted as parameter changes due to evolution, and the probability of 
exiting after t steps represents the probability that the network is no longer capable of correctly performing 
its function (for instance, when a lethal mutation occurs). Studying the first exit problem is the natural thing 
to do in certain evolutionary models. Suppose we consider a fitness landscape on the parameter space where 
the functioning regions have a fixed high fitness and every other region has zero fitness. If we consider a 
space of alleles to be nearly continuous and model the effect of mutation as diffusion in this space, as is 
often done in the adaptive dynamics literature 1241 . we find that we need to compute the mutation load, 
namely the rate of death from exiting the high fitness region. This idea was previously used in the context 
of transcriptional networks [4 |. 

To explore the shape of the regions Gc.i to Gamo, Algorithm I uses a random walk in the parameter space, 
and checks "exit times" as well as the "failed parameters". 

Algorithm I 

Pick a positive number oq to be the constant magnitude of the random walk step. 
Repeat points 1-4 (run q), Q times. 

1. Step 0: generate a point = {p\, . . . , p^)' at random in the parameter region G-^; 

2. Step A; — 1/2, A; > 1: generate a random perturbation dp^ G [— «0) «o]™> such that \dp^\ = aoQ 

3. Step k, k > 1: check if p^ = p'^^^ + dp^ is still in G-^; 

4. Check. The random walk exits the parameter region at time t if p^ G for A; < f but ^ G^. 
Let p*-^ , . . . , ^ be parameters that fail to satisfy the hierarchy of conditions which defines G^. 
Update the exit times vector: exit (q)=t. 

Update the failed parameters vector: f ailpar {q) =[ji, . . . , jj]. 



'This step corresponds to generating a random point from a uniform distribution over the hypersphere in n dimensions, which 
can be achieved by the Box and Muller transformation 1 25 1. Briefly, for i = 1, . . . ,n pick Zi randomly from a gaussian distribution 
of mean zero and variance one. Then normalize to obtain z — (zi, . . . , Zn)/ \/ z\ + • ■ ■ + z'^. 
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To interpret the numerical results obtained with Algorithm I, define the probability that a mutation takes 
place in the first t steps by: 

P.Ut) = ^ card(/t), It = {qen: exit {q) < t}, 

where card(-) denotes the cardinality of a set. Algorithm I was applied to each component of the feasible 
parameter space of the segment polarity network, with oq = 1 x 10^"^ and Q = 4000. Two striking facts are 
revealed. First, with a significant probability, fluctuations in the parameters will drive the system from the 
operating regions Gci, Gcu, Gem, or Gciv, to the region Gauio and, conversely, switching was also observed 
from Gauio to the other four (see the switching column in Table |4]l. Recalling the difference between Gamo 
and the other four components, this means that, in a significant number cases, the network responds to 
perturbations by switching to an alternative biological pathway, rather than break down. A second fact is 
that only a very small number of parameters (six out of 39, namely Kcn^,,, Kwg™, Kci,„c, Kpjca, ^ciug> ^-cNwg) are 
responsible for above 90% of network failures or mutations. The percentage of cases where each of these 
parameters failed is shown in Table |4] 



Table 4: Fragile parameters. 



G^ 


% switching 


% failed 


% failed 


% failed 


% failed 


% failed 


% failed 


















Gc,i 


1.21%(toGA„J 


5.4% 


6.3% 


20.0% 


4.4 


50.8% 


5.6% 


Gc,U 


0.88% (to GA„,e) 


4.7% 


5.6% 


17.2% 


4.3 


40.0% 


18.7% 


Gc,lU 


1.50% (to Ga„,o) 


6.2% 


5.2% 


30.8% 


4.1 


35.0% 


11.5% 


Gc.lW 


1.16% (to Ga„.o) 


4.9% 


5.3% 


16.8% 


5.0 


40.1% 


18.9% 


GauIo 


0.02% (to Gci-Gciv) 


8.5% 


5.7% 


21.2% 


4.8 


31.4% 


18.8% 



Calculating the distribution function i^^„, shows that the probability of mutation increases very rapidly 
for small times, in all five components (see Fig. [6]) - this indicates a low robustness of the network, because 
it is very likely that a very small number of fluctuations leads out of the feasible parameter space. To 
compare the results for the five components, we computed some quantities of interest. A possible indicator 
of robustness is T1/2, defined as the time for which there is a 50% chance that the system has already suffered 
a mutation. Low T1/2 indicates a system which has a low robustness to perturbations. Another indicator 
is -Pn,i„(10), which gives the probability that the system has been disrupted after only 10 perturbation steps. 
Similarly, the values P„,^t{WO), Pn,„,(1000), and P„,„,(10000) are also shown for comparison. The computed 
values are summarized in Tablepl Comparison of the values for and Pn,u,(10^) {d = 1, . . . , 4) in the five 



Table 5: Indicators of robustness. 



G-y 


Tl/2 






P„.,(1000) 


P„.,(10000) 


Gc,i 


18 


0.40 


0.74 


0.94 


0.99 


Gc,U 


23 


0.37 


0.72 


0.94 


0.99 


Gc,m 


19 


0.40 


0.72 


0.93 


0.99 


Gc.IV 


23 


0.37 


0.73 


0.94 


0.99 


GauIo 


309 


0.30 


0.41 


0.61 


0.88 



regions, shows clearly that the components Gc,i to Gc.iv result in a less robust network, while Gauio exhibits a 
much higher level of robustness. Furthermore, there is a non-negligible probability (« 1%) that the network 
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switches from the other components to Gau,o instead of breaking down, thus contributing to the robustness 
level of this component. Distinguishing the levels of robustness among the four regions Ga to Gc,iv is not 
straightforward, since the indicators P„,i„(10'^) and T1/2 are very similar. 

As noted above, only six out of 39 parameters are responsible for over 90% of failures. Curiously, two 
of these parameters satisfies constraints which are in fact independent of the regions (kci^k and Kptcci> see 
Table [7]). So, the numerical results clearly show that the parameter space is very narrow in the directions 
defined by Kci,„c and Kptcci- The other critical directions are defined by Kq,,,, and KcN„g, and Kcn™ and K^oe,,, 
which satisfy different conditions in each of the five parameter regions (Table [8]). Together with common 
parameter Hcipu-, the parameters that regulate activation and inhibition of wingless by Cubitus proteins (kuh^ 
and Kctiwg) are the most critical. 

The main conclusion from Algorithm I clearly follows the preliminary estimates of the relative volumes 
(compare Tables [T] and |5] both concluding that Gauio is more robust than the other four components). But 
the geometry analysis reveals three new fundamental results: (i) the system increases its robustness to envi- 
ronment perturbations by switching to an alternative biological pathway. The switching event may be from 
a "small" to a "large" region but also, more remarkably, from a "large" to a "small"; (ii) the lack of robust- 
ness is due not only to small sized regions, but in part to critical parameters {kqipk and Kptcci)> which define 
directions along which the parameter space is globally very narrow; (iii) the volume alone is not a reliable 
measure of robustness, since volume (Tablejljl and the indicator T112 provide different robustness classifica- 
tions for components Gc,i to Gciv For instance, the volume of Gen is apparently the smallest (an indicator 
of low robustness), but is the largest (an indicator of high robustness), suggesting that the shape of the 
region does plays an important role. In contrast, the numbers P,^„i{\{)'^) are very similar, suggesting that 
robustness levels of Gc,i to Gc w are in fact very similar. However, it should be noted that neither volume nor 
Tij2 provide conclusive information on the relative levels of robustness of Gcj to Gc,i\- In particular, note 
that Ti 12 depends on the magnitude of the random walk step - other numerical experiments were performed 
with different oq values (not shown), and the comparison results are unchanged. 



5 Discussion and conclusions 

Analysis of the feasible parameter set, by estimating its volume, identifying connected components, and its 
geometric properties aie valuable tools for establishing and quantifying robustness in regulatory networks. 
The concept of robustness, in the sense that the system's regulatory functions should operate correctly under 
a variety of situations, is closely related to the parameter space and the effect of parameter perturbations. In 
this context, our analysis suggests that the segment polarity network is vulnerable to perturbations in its pa- 
rameters. Indeed, the first striking result from our analysis is that the feasible parameter space is composed 
of five disconnected components. An implication of this topological characterization is a diminished ca- 
pacity of the network to respond well to environmental perturbations. Random fluctuations will often drive 
the system to a set of parameters outside any feasible region, and thus lead to a break down of the network 
or a different phenotype. Indeed, as the results of Algorithm I show, sucessive random perturbations to the 
parameters will drive the system out of the feasible parameter set, with a large probability. For instance, if 
parameters are randomly perturbed for up to 10 times, each of magnitude 1 x 10~^ in any direction, there 
is a 30% probability that the system will fail to operate correctly (see Table |5} column P„,„,(10)). On the 
other hand, it is possible that a series of fluctuations in the environment may drive the system to adopt an 
alternative biochemical pathway, and thus "jump" from one feasible component to another (with probability 
1%, seeTable|4]). 

As the group of most fragile parameters suggest, the Qv!o\y\x%-wingless interactions are at the basis of 
the appearance of disconnected regions of parameters. Dis-connectivity in the space of parameters can be 
traced in large part to an incompatibility of Cubitus repression functions in the second cell: CN2 should 
be present to repress engrailed expression, but should be absent to enhance CI2 activation of wingless. To 
increase the network's robustness to environmental fluctuations, the segment polarity model should account 
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for engrailed regulation by other factor than Cubitus. One possibility is to include regulation by pair-rule 
gene products, such as Sloppy paired, as explored both in [9| and [8|. An external factor, again possibly 
from the pair-rule genes, will also play a major role in establishing asymmetry in the cubitus levels (Ui). 
These contribute to a larger admissible parameter space, and together with an improved engrailed regulation, 
will greatly enhance robustness of the segment polarity network in maintaining its pattern. An extension of 
the current analysis including the regulation by Sloppy paired is currently in preparation by A. Dayarian at 
one of our labs ||26]| . 

Both the volume estimates and the probability of failure or mutation (P,^^) in each component indicate 
that Gauco is the most robust parameter region, while Gci, Gc.n, Gc.m and Gc.iv are less robust regions, all 
at the same level. However, volume is not a reliable indicator of robustness by itself, and fails to predict 
alternative robustness mechanisms. Additional knowledge on the network mechanisms has been gained 
with the geometry analysis. A noteworthy fact is the non-negligible probability (1%) that fluctuations in 
the parameters in regions Gc.i to Gc.iv result in a switch to the region Gp,^to, and remarkably (but with lower 
probability 0.02%) also from Gauio to the others. Of the five disconnected components, Gc.i, Gen. Gc,m, and 
Gciv correspond to the pathway where wingless is regulated by Cubitus interruptus proteins, while Gamo 
corresponds to the pathway where wingless is regulated by its own protein levels. Thus it is more likely 
that wild type expression in the segment polarity network is achieved through the Wingless auto- activation 
pathway. In the absence of the auto- activation link, von Dassow et. al. failed to observe any feasible 
parameter set in their numerical experiments. However, as soon as the auto-activation pathway was added 
(the second "missing link" in the model Q), immediately a significant percentage of feasible parameter sets 
were observed. This is not surprising, as elucidated by our analysis: while wingless auto-activation is not 
strictly necessary to establishing the segment polarity genes pattern, it does greatly increase the probably 
that the pattern is achieved (Gauio has a much larger volume, by a factor at least 40, and also exhibits higher 
robustness indices). 

Another fundamental conclusion from the geometry analysis is the existence of six (out of 39) critical 
parameters which are responsible for 90% of the network failures due to parameter fluctuations. Moreover, 
the intervals for two of these parameters (kci;« and Kpjca, Table|7]l are independent of parameter space com- 
ponents. The feasible parameter set is thus globally restricted by these parameters, which define "narrow" 
directions (see Fig.[T](b) ). 

Robustness of a regulatory module should not be measured simply as a function of the volume of its 
admissible parameter space. The geometry (for instance, convexity or existence of sharp points) and topol- 
ogy (connectedness) of the parameter space play fundamental roles in measuring robustness. The analysis 
developed in this paper can be applied to other systems and regulatory networks, to systematically charac- 
terize and explore the admissible space of parameters, its topology and geometry. These provide reliable 
information on how the network's interactions contribute to its robustness or fragility, and serve as measures 
to classify robust regulatory modules. 
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A Notation 



The original model can be found in ||3]|l3l. In order to make our work more clear, we include the notation 
as well as the original equations below. Without loss of generality (the geometry remains unchanged), 
each cell is assumed to have four faces (Fig. [2]), rather than six as in the original model [3|. The model 
reproduces a parasegment of four cells and uses repetition of this group of four cells to reproduce the 
embryo's anterior/posterior axis (A/P axis in Fig.|2]l, and the circular ventral/dorsal axis (V/D axis in Fig. [2]). 
Because intercellular diffusion is only considered along the A/P axis (left/right), and because cells repeat in 
the orthogonal V/D direction (up/down), it is indeed equivalent to consider symmetric four-sided or six-sided 
hexagonal cells. 

I Four— cell parasegment ■ 



V/D 



wg 



wg 



A/P 



wg 



en 



en 



en 



wg 



en 



4- wg 1 4- en 1 



4- wg 1 4- en 1 



Figure 2: Four cells in a parasegment, with periodic boundary conditions in both dimensions. Each cell has 
four membranes. The relative values of Wingless in each cell (EWGj) are shown. 

A saturation function, and its horizontal reflexion, are introduced: 



1 — 4'{X^ K, u). 



(I){X, K, V) 
IpiX, K, V) 

The subscripted variables are as follows: 

Xi = concentration of species X on cell i (when homogeneous throughout the cell ), 



X 

n{i,j) 

-^n{i,j),j+3 
X: 



concentration of species X on cell i, at face j, 
threshold for activation of species Y, induced by species X, 
index of neighbor to cell i, at face j, 
concentration of species X on cell face apposite to i, j, 

6 

Xij = total concentration of species X on cell i, 

6 



total concentration of species X presented to cell i by its neighbors. 
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B Original equations 

From (Si [131, the model equations are: 

((^(EWGiV'(CNi, Kctien, I^CN™), Ky^Oe,,, '^WGen) " etli) 

-{erii - ENj) 

1 



derii 


1 


dt 


Hen 




1 


dt 


Hen 


dwgi _ 


1 


dt 





dt 

dEWG 



dptci 
dt 

dPTCjj 

dt 



{Wgi - IWGj + Tend^ifiwoEWGij - iJwoJ^exoIWGi) 

re,„IWGi - r,„d„EWGjj- + r-Af(EWG„(jj)j+3 - EWGj^j) 

EWG, 



1 
6 



+rLM(EWGij_i +EWGij+i - 2EWGij) - — ^ 



1 

1 /I 



) 'I'CNpfc) i^CN/irc)) ^Clptci 

ptCi — VTCi j — ft;pTCHH-f^PTc[HH]oHH„(j ,,_|_3PTCjj 



dcii 


1 ^ 


~dt 


Hci 


dCli 


1 


dt 


Hci 




1 


dt 


Hci 


dhhi 


1 


dt 


Hi,i, 




1 


dt 


Huh 



HpYc \ 6 

+rLA/PTc(PTCij_i + PTCij+i - 2PTCij) 

(c/j — Clj — i?ciCciCIj(/)(PTCj J, KpTCCIi '^PTCCl)) 
(i?ClCciCIj(^(PTC},j, KpTCCI) J^PTCCl) ~ CNj) 

(i;^(ENj'i/;(CNj, Kcnwd ^cnaa)) ^en/i/d z^enwi) ~ /j^j) 

/j/jj — HHjj — Kp-rcHH-ffHH[PTC]oPTC„(j j) j_|_3HHj.j,- 
+?'LMHH(HHjj_i +HHjj+i — 2HHjj) 
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Simplified model, for large u 



fern = {di^y^Gil^WGenWil^CNen " CNj) - CTj) 

/sNi = ^{erii-ENi) 

-"EN 



1 / QciH.g^(CIt - /tciH.g)6'(/tcNH.g - CNj) + awGM.g^(IWGi - KwgwJ 



/iWGi = 7^ (Wgj - IWGj + rendo-fflWoEWGi^T - -f/wG?^exoIWGi) 
-"WG 

/sWGij = l^exoIWGi - TendoEWGi J + rM(EWG„(j j)_j+3 - EWG^ j'^ 

EWG 

+rLM(EWG,,,_i + EWGi j+i - 2EWGi j) - ''^ 



WG 



fp,ci = {Tie{Cli - Kap,MKa,p,c - CNi) - ptCi) 

+rLMPTc(PTQj_i + PTCij+i - 2PTCij) 
fcti = ^ {0{Bi - Kj,„)9{Kw,i - ENi) - cii) 

fcH = (di - Cli - ifciCciCIie(PTQ,T - KPTCCI)) 

/cNj = (-f^ClC'ciCIi^(PTCj^T ~ /^PTCCl) ~ CNj) 

A^i = (6'(ENi - Kenm)^(kcn«. - CNi) - hhi) 

fwiij = ^^/J^i — HHjj — KpTCHH-f^^fffl[PTC]oPTC„(j j+aHHj^ 

+?'LMHH(HHjj_i + HHij+1 — 2HHij) 
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D Steady state pattern 



Solving equations (31l-(42i at steady state (/ = 0), and simplifying where possible, yields the algebraic 
expressions: 

eui = ^(EWG, - KwG.„) ^(kcn™ - CNi) (43) 

ENi = em (44) 

^ _ Qci.rg^(CIi - Kci,,g)g(ft;cN,.g - CNj) + awGH.g^(IWGi - KwowJ ^^^^ 

1 + aci„,g^(CIj - Kci,vs)6'(KcNvrs " CNj) + OwG.vg^' (iWGj - KwGws) 

IWG. = ^-^ore. EWG,^+ ^ wg. (46) 

1 -|- IJ-Ywa' exo J- ~r -fliwG' exo 

MEWG = v?§ (47) 

4 1+ -f/lWG?'exo 

= r,e(CIi-Kc,„.)e(ACcNp,c-CN,) (48) 

PTCjj = ^P^Cj — KpTCHH-f^PTc[HH]oHH„(j ,,_(_3PTCjj 

+rLMPTci^PTc(PTC,,,_i + PTC,,,+i - 2PTC,j) (49) 

cii = [/, 0(ken„ - ENi) (50) 

CI, = ^(^ENc-ENQ ^^^^ 

1 + HqiCqi ^(PTCj.T '^PTCCl) 

CNi = Ui ^ ^"f " g(AtEN., - EN,) g(PTC,,T - KpTcci) (52) 

hhi = ^(ENi - Ken,,,,) ^(kcn«, - CNi) (53) 



1 



-'^ij ~ ^/^^j — ^PTCHH-f^HH[PTC]oPTC„(jj)j_|_3HH, 



+rLMHHi?HH(HHi,,_i + HHi,,+i - 2HHij) (54) 

EWG is a vector in M^^ with components: 

EWG = ( EWGi,i,EWGi,2,EWGi,3,EWGi,4,EWG2,i,EWG2,2,EWG2,3,EWG2,4, 
EWGs,!, EWG3,2, EWG3,3, EWG3,4, EWG4,i, EWG4,2, EWG4,3, EWG4,4)' 

wg is also a vector in M^^, given by the following Kronecker tensor product 

wg = (wgi,Wg2,Wg3,Wg4)' Xfcron (1,1,1,1)' 

= {Wg^, Wg^ ,Wg^,Wg^, Wg2 , Wg2 , Wg2 , Wg2 , Wg^ , Wgr^ , Wgr^ , Wgr^ , Wg^ , Wg^ , Wg^ ,WgJ. 
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Putting together the 16 equations (35 1, and substituting IWGj by its steady state expression (34i, it is not 

; M}^ is composed of various 4x4 blocks, as follows: 



difficult to see that the matrix M G 



M 



where 



E 



with 



/ E 


-^24 





Fa2\ 




E 


-P24 








F42 


E 


F24 


\F2A 





F42 





f-d 


rLM 


rM 


rLM\ 


rLM 


-d 


rLM 







TLM 


-d 


rLM 


\rLM 





rLM 


-dj 



+ h 



A 


1 


1 




1 


1 


1 


1 


1 


1 


1 


1 


V 


1 


1 


1/ 



d- 
h 



-1 



IWG ' endo + rM + 2rLM, 

^endo 



(55) 



IWG' cxo 



^24 



/o 








o\ 











rM 














^0 








/ 



M2 



F' 



24 











0^ 


























Vo 


rM 








Note that the steady state equations for EN, IWG, EWG and PTC are algebraic, and in fact exact solutions 
can be computed from the steady state values of wg and ptc. These are discussed in more detail in the 
Appendices |E] and |F] 



Remark: The parameters are as in [3], except Ti and Ui, which represent the maximal values of ptc and 
ci (respectively), in each cell. These take values in the interval [0, 1] and generalize the possible ON values 
of ptc and ci. 



Note that, in the simplification from ( [24| ) to p6| ), we have generalized the equation and added distinct 
maximal levels of expression in each cell, given by Tj (i = 1, . . . , 4). This allows a more accurate represen- 
tation of the experimental, which shows that patched is strongly expressed in every second and fourth cells, 
weakly expressed in every first cell, and not expressed in every third cell (see [3| for more discussion). Thus 
we will consider the case: Ti <T2 = T4. 

A similar generalization was made to deal with the activation of cubitus interruptus. In von Dassow 
et. al. model, this is due to some external parameters Bi (not governed by a dynamical equation), with a 
corresponding activity threshold Kb«- However, for more generality, and to allow distinct maximal levels 
of expression in each cell, we have replaced each of the terms 9{Bi — k^,) in (38 1 by a parameter Ui, 
i = 1, . . . ,4 (50l. Furthermore, in characterizing the set of feasible parameters, it will become clear that 
allowing distinct Ui enlarges the space of possible parameters, by introducing the four regions Gq i to Gciv 



E Analytically solving Wingless levels 



The steady states of Wingless proteins (47i and (46 1 are given directly by algebraic equations, depending 
only on wingless mRNA {wg2) and diffusion parameters for intracellular (membrane-to-membrane) and 



intercellular communication. Consider equation (47 1: it is easy to see that M is in fact always invertible (if 
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all parameters are positive). First note that the matrix is diagonally dominant, by adding up the entries in 
any column: 

- {H~l + r,„j„ + rM + 2rL]\i) + 2rLM + tm + Ah = -H~\ - r,„j„ ^ 

which is always a negative quantity. By Gersgorin's Theorem, all eigenvalues of M are contained in the disk 
centered at —d + h with radius 2rLM + fM + 3/i, therefore all have negative real parts. Thus, the matrix M 
is symmetric and negative definite, and since the right-hand-side vector is also non-positive, all solutions are 
real and positive, whatever the choice of parameters. As a fact, note that the vector 1 = (1, 1, . . . , 1)' G M^^ 
is an eigenvector of M, corresponding to the eigenvalue Ai = —H^\ — r^ndo u^^^j^j. — ■ 

Proof of Theorem|2] Assume that wg = (0, w, 0, 0), for any positive constant w. From the symmetry of 
the matrix equation ([47]), several facts can be deduced, which lead to the main result ([T3]). 



Fact E.l. For alH = 1, 2, 3, 4 it holds that 

EWGi,i = EWG,,3. 
Proof. This is easy to see from the respective equations: 

{-d + /i)EWGi,i + {tm + /i)EWG,,2 + {rM + /i)EWGi,3 + (^m + /i)EWGi,4 = J^'^Si 

(rM + /i)EWGi,i + (rM + /i)EWG,,2 + {-d + /i)EWGi,3 + {tm + /i)EWGi,4 = J^—'^Si 

which can be rearranged to 

-{d + rM)EWGi,i + (rM + /i)(EWG,,2 + EWGi,4) + {ru + /i)(EWGi,3 + EWGi,i) = ^>vft 

-{d + rM)EWGi,3 + {tm + /i)(EWG,,2 + EWGi,4) + {ru + /i)(EWGi,3 + EWGi,i) = %—w8^ ■ (56) 

?"endo-"IWG 

Subtracting these two equations yields the desired result. I 
Fact E.2. It holds that 

EWG2,2 = EWG2,4, EWG4,2 = EWG4,4, EWGi,2 = EWG3,4, EWGi,4 = EWG3,2. 
Proof. Exchanging the indexes: 

2,2^2,4 4,2^4,4 1,2^3,4 1,4^3,2 
it is easy to see that the system remains unchanged (see also Fig.|2]). I 



The equality part in ( 1 3 1 is now clear: 
Fact E.3. EWGi = EWG3. 

Proof. We first show that EWGi^i = EWG3^3. Writing equation ([56]) for i = 1 and i = 3: 



-{d + rM)EWGi,i + (rM + /i)(EWGi,2 + EWGi,4) + (rM + /i)(EWGi,3 + EWGi,i) = 
-{d + rM)EWG3,3 + {tm + /i)(EWG3,2 + EWG3,4) + {tm + /i)(EWG3,3 + EWG3,i) = 
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Using Fact E.l one has EWGi^i = EWGi^a and EWGs^i = EWG3 3, and then using Fact E.2 obtains: 

-{d + rM - 2rM - 2/i)EWGi,i + {vm + /i)(EWG3,4 + EWG3,2) = 
-{d + rM - 2rM - 2/i)EWG3,3 + {tm + /i)(EWG3,2 + EWG3,4) = 0. 

Subtracting these two equations shows that EWGi^i = EWG3,3. Now recalling the notation for Xj from 
Appendix [A| 

EWGi = EWGi,i + EWG2,4 + EWGi,3 + EWG4,2 
EWG3 = EWG3 1+EWG4 4 + EWG3 3+EWG2 2. 



Using EWGi,i = EWG3,3, Fact E.l and Fact E.2 obtains: 

EWGi = EWG3,i + EWG2,2 + EWG3,3 + EWG4,4 = EWG3. 
as we wanted to prove. 

To show the other inequalities, note first that the 16 variables EWGjj are thus reduced to only seven: 

3,3 



-£■1,1 


= EWGi,i 


= EWGi, 3 = EWG3,i = EWG 


-El, 2 


= EWGi, 2 


= EWG3,4 


-£-1,4 


= EWGi, 4 


= EWG3,2 


-£^2,1 


= EWG2,i 


= EWG2,3 


E2,2 


= EWG2,2 


= EWG2,4 


-£■4,1 


= EWG4,i 


= EWG4,3 


-^4,2 


= EWG4,2 


= EWG4,4 



and satisfy the equations: 



{d-rM-2h)Ei,i + {rLM + h)Ei,2 + irLM + h)Ei^4 = (57) 

2{rLM + h)Ei,i-{d-h)Ei,2 + hEi,i + rME2,2 = (58) 

2{rLM + h)Ei,i + hEi^2-{d-h)Ei,4 + rME4,2 = (59) 

-id-rM-2h)E2,i+2{rLM + h)E2,2 = %—wg2 (60) 

h 

2{rLM + h)E2,i-{d-2h)E2,2 + rMEi^2 = 73—^^2 (61) 

-{d-rM-2h)E4,i+2{rLM + h)E4,2 = (62) 

2{rLM + h)E4,i-{d-2h)E4,2 + rMEi,4 = 0. (63) 



To simplify notation, set: 



A = d-VM -2h, B = 2{rLM + h), w= ^- — wg2, 



and note that A > B > 0. 
Fact E.4. The following hold: 

(a) £'4,1 < £^4,2 < -E'1,4 < -E'1,2 < -E'2,2; 

(b) Eai < E,i < Eix, 
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(c) £'1^2 + -E'1,4 < -£'2,2 + -E'4,2 



Proof. To prove part (a), from eqs. ( 62 1, ( 63 1 it holds that 



F -^F ■ F 

^4,l-^i?4,2, ^4,2 - ^2 _ 52 + 



E 



1,4 



Because ^4 > i? > 0, it is clear that i?4^i < -£4^2 < -£^1,4- From eqs. (60 1, (61 1 it holds that 
B 



w 



Then eqs. (58 1, (59 1 can be written in the form 

TmA 



d - VM 

d-TM 



A'^ -B^ + tmA 

ruA 
A^ - B"^ + TmA 



El; 



BEi^i + h{Ei^2 + £1,4) + VM 



A + B 



A^ -B"^ + vmA 



w 



-El, 4 = BEii + h{Ei^2 + -El, 4) 



which implies that E'l ^ < £'1,2 (it is easy to see that the factor multiplying both £1,2 and £1,4 is positive, 
since d > tm)- 



We still need to prove the last inequality in (a), but we can now prove (b). From eq. (57 1 

-^1,1 = 2 ^4 ^^^'2 ~^ -^1,4) < -£^1,2 
using (a) and because B < A. This proves the second inequality in (b). To prove (c), substitute this £1,1 

(-£1,2 + -Ei,4) > -£1,2 + El,4- 



expression into the sum of eqs. (58 1, (59i: 

A^-B^ + tmA 



£2,2 + £4,2 



TmA 



The last part of (a) now follows from (c) together with £4,2 < El, 4, which implies £1,2 < -£2,2- 
Finally, the first part of (b) is easy to see from: 

1 B B 1 B 



To prove the first inequality of Theorem [2] is now straighforward. 
Fact E.5. EWG4 < EWGi 



Proof. Recall the notation for EWG^ and use Fact E.4 
EWGi - EWG4 



2Ei^i + £2,2 + £4,2 — 2^4,1 — 2E'i,4 
2(El,l — -£4,1) + (-E'2,2 + -^4,2 — -^1,2 — -E'1,4) + (-El, 2 — -^1,4) > 0. 



The next result finishes the proof of Theorem |2] 
Fact E.6. EWG3 < EWG2 
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Proof. Consider: 

EWG2 - EWGi = 2^2,1 + 2^1,2 - 2^1,1 - -^2,2 - Ea,2 

= 2(Si^2 — Ei^i) + (-52,1 — -£'2,2) + (-E'2,1 — -£^4,2), 

which is positive if £^2,1 > £^2,2- We will show that this is indeed the case. It will be useful to see that 

A + B 1 A^-B^ + vmA -A^ + B^ + {d- rM)A _ 

i'2 ^^^d(yi2 _ 52 + ^^^^) -rl.A'^^ 2''^'d{A^ - 52 + rMA) - t\A {A - B){A^ - ^2 + ^ruA)""- 



Now consider 



A-B tmA A-B A^B _1 

A A'^-B'^ + rMA ^'^ A~A'^-B'^ + rMA'^^A'' 



The last two terms can be combined into 

rM 



w, 



-B'^ + tmA 
and the two terms due to i?i 2 can be simplified to: 

1 {A-B){A + B) 

^A^^2 _ 52 + ^^^^ j_^^2 _B2 + rMA) - tmA"" 

and 

TM 1 -A^ + B'^ + {d- rM)A 



2 A2 - 52 + 2rMA jL^a^ _ 52 + _ ^^jA 

Factoring out rMw/{:^{A'^ — B"^ + vmA) — vmA), one obtains 
.(^2 _ ^2 ^ ^^^^^^ _ rMA){E2,i - £2,2) = 



w. 



■-rM 

1 d 



TMW rM 



' ( \a--B- + rMA) - rMA - (A^ - B^^^ ' "^^ + + (d - r.M 



A^-B^ + rMA\rM J 2 A^-B^ + 2rMA 

which can be further simplified to 

4-1) - ^ 1 A^-B^ ^ {d-rM)A _ 1 {d-rM)A ^ ^ 



A^-B'^ + rMA 2A^-B'^ + 2rMA A^-B'^ + rMA 2 - B"^ + 2rMA 

because the first two terms are clearly positive, and the last two terms add up to a positive number. This 
shows that £'2,1 > £2,2, as we wanted to prove. I 



F Analytically solving PTC and HH levels 

In this section, we prove uniqueness of solutions for PTC and HH in the conditions of Theorem [T] The 
steady state levels of Patched and Hedgehog proteins are given by a system of nonlinear equations ([49]) 



and (54 1. These equations can be solved explicitly and uniquely in the case ptC2 = ptc^ = T2, which is true 
is the steady state output is in 3^*^. To simplify notation, we use 

rp = tlmptc, rn = riMm, i^H = Kptchh[HH]o, up = Kptchh[PTC]o, 
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and define 



1 , 1 

"P = 77 V2rp, dH = jj: \-2rH. 

-"PTC -"HH 

We introduce furtiier notation: 

_ 2rj,dp _ 1 1 2rp{rp + dp)_ 1 dp{2rp + dp) 

dl- 2rl ' ~ ^ 4/7pTc (fp - 2rl ~ 4H„c - 2rl ' 

Lemma F.l. Let x G W be sucii tliat h{x) G 3^*^. Tiien, tiie solution for HH is: 

HHj 1 = HHj 2 = HHj 3 = HHj 4 = 0, i = 1, 2, 4, 
HH3,2 = HH3,4 = Root+ , 

HH3,i = HH3,3 = :^(7^ + rHHH3,2 + rHHH3,4), 
where Root+ is the positive root of the quadratic equation: 

knidjj - 4rl)X^ + (^{dp - Pp){d]j - Arl) - knidn + "^^h)^ + dnkp-fpptc^ X 

-{dp-(3p){dH + 2rH)^ = 0. 

4-nHH 

And the solution for PTC is: 

PTC3,1 = PTC3,2 = PTC3,3 = PTC3,4 = 0, 
PTC2,2 = PTC4,4 - 



dp - (5p + kHYai^,A 

PTC2,i = PTC2,3 = PTC4,i = PTC4,3 = \ 2 irpdpVTC2,2 + tt^C^p + rp)T^ , 

dp - 2rf, V 4iipTc / 

PTC2,4 = PTC4,2 = -^(7^+2rpPTC2,i). 

dp 4 iipTc 

Proof. Let x G W and h{x) be a vector in y"^^, defined by ([t]). Because hedgehog is not expressed in 
cells 1, 2 and 4, note that for i = 1,2,4 

4 4 

HH,,T = = hh, - Kp{- ■■)+rHY, (HHij^i + HH,j„i - 2HH,j) 

= -Kp(---) 

since /z/j = (0, 0, 1,0), and the sum that multiplies th cancels out. The terms in Kp(- • • ) are all nonnegative, 
and therefore they can only be zero. We conclude that: 

HHi,i = HHi,2 = HHi,3 = HHi,4 = 0, i = 1, 2, 4. 

A similar argument shows that ptc^ = implies: 

PTC3,1 = PTC3,2 = PTC3,3 = PTC3,4 = 0. 
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Therefore, the only nonhnear terms appear in the equations for PTC2,2 and PTC4,4: 

dpPTC2,2 - rpPTC2,i - rpPTC2,3 + KffPTC2,2HH3,4 = -r^ptC2 

4-npTc 

dpPTC4,4 - rpPTC4,i - rpPTC4,3 + At/fPTC4,4HH3,2 = -rn—ptCi ■ 

4-npTc 

Moreover, symmetry of the system shows that PTC2,i = PTC2,3 and PTC4 1 = PTC4,3, because each pair 
satisfies exactly the same equation: 

a!pPTC2,i - rpPTC2,2 - rpPTC2,4 = j^PtC2 (64) 
dpPTC4,3 - rpPTC4,4 - rpPTC4,2 = j^Ptc^. (65) 

We then have: 

PTC2,4 = ^(^^^^'^2 + 2rpPTC2,i) 
PTC4,2 = ^(i^^^^4 + 2rpPTC4,i). 

Solving for PTC2,i as a function of PTC2,2> and for PTC4,i as a function of PTC4^4: 
PTC2, = 4^(^P^pPTC2,2 + ^(c^P + rp>.C2) 
PTC4,i = ,2 \ 2 ( r-pdpPTC4,4 + -7^{dp + rp)pto4 ) . 

Thus we get equations depending only on PTC2,2 and HH3,4, and on PTC4,4 and HH3,2: 

dpPTC2,2 - ,2 o 2 ( ?^pc^pPTC2,2 + -r^{dp + rp)ptcA + kjyPTC2,2HH3,4 = j^PtC2 (66) 
dp — ZVp \ 4iipxc / 4jipxc 



2rp / 1 \ 1 

,2 o 2 rpcipPTC4,4 + -rE—{dp + rp)ptc^ + Ki^PTC4,4HH3,2 = — 

Up — 2,T r} \ 4x2 PTC / 4x1 p' 



On the other hand, since PTC3J = for all j, it follows that: 

HH3,i = HH3,3 = ^{j^hh3 + ri/HH3,2 + rj^HH3,4), 
and substituting into the HH3 4 and HH3 2 equations: 

Tf^ 1 1 

d//HH3,4 - 2-^(77^^3 + rHHH3,2 + rj^HH3,4) - KpPTC2,2HH3,4 = j^^hh (68) 

rpf 1 1 
rfffHH3,2 - 2 -/^ (77^/1^3 + r/^HH3,2 + ri/HH3,4) - ACpPTC4,4HH3,2 = jr—hhs . (69) 

"if 4iiHH 4iiHH 
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The last four equations may be solved for the four variables PTC2,2» PTC4^4 HH3 2 and HH3 4, and the 
remaining PTC, HH will then follow. Recalling the notation introduced above, one can write 



dp- I3p + fcffHH3,4 dp- I3p + fcHHH3,2 



This leads to 



d//HH3,4 - 2^(-^— /j/j3 + rHHH3,2 + r//HH3,4) - Kp-fpptC2-, J^^^f = T^^^s 

dn 4i?HH dp- I3p + fc/fHH3,4 AHna 

ru , I ^ HH32 1 

a//HH3,2 - 2:^(773— ««3 + ?-hHH3,2 + r//HH3,4) - Kpjpptc^- — , ) = jrj—hh ■ 

dn 4i?HH dp - (3p + A:/^HH3,2 4iJHH 

From the symmetry of these equations, it is easy to see that 

ptC2 =ptCi =^ = HH3,2- 

and thus have the following equation for HH3 4 = HH3 2 = X (after some simple algebra steps): 

knidl - 4rjj)X^ + (^{dp - Pp){dl - Arl) - knidn + '^^h)^ + dukp^pptc^ X 

-{dp-/3p){dH + 2rH)^ = 0. (71) 

4-nHH 

We next show that only one of the two roots of this second order polynomial is positive and hence the unique 
solution to HH3 2, HH3 4. Let the polynomial be of the form 02^^ + ciX + cq = 0. The term inside the 
square root will be cf — 4coC2 where: 



-4coC2 = 4kH{dij - 4rij){dp - Pp){dH + 2rjy)- 

4-nH 

The factor d'j^ — Ar'jj is positive, by definition of dn- The factor 

irl^-^^/^2 ,^2 



dp- (3p = dp- ,2 _n 2 = '^P^dp - ^^P) 



is also positive, again by definition of dp. This means that c\ — 4coC2 > cf , so whatever the sign of ci, 
— ci — ■\/c\— 4coC2 < 0, which leaves us with: 

nH3 2 — HM3 4 — 

^C2 



(the coefficients are as in (71 1). 



El Asymmetry in patched ON levels 

The assumption T2 = T4 is now relaxed, and the more general case is analyzed. The main question is how 
Patched asymmetry influences the space of parameters, G, and whether the five components can become 
connected. In other words, does the more general case assumption T2 / T4 leads to a increasing network 
robustness. It will be seen that this is actually not true. The presence of CN in the first cell is still necessary 
(because Wingless protein expression is not affected by ptc levels), but expression of CN in the second and 
fourth cells may now be different. While it is now difficult to explicitly solve the nonhnear equations for 
PTCj and HHj, it can still be shown that ptC2 < ptc^ implies PTC2 < PTC4. 
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FactF.l. {ptC2 -/7te4)(PTC2 -PTC4) > 0. 

Proof. To see this assume that ptC2 > ptc^ (the opposite case follows a similar argument). From the 
discussion above, the Hedgehog values must satisfy 

HH3 4 HH3 2 

dHHH3 4 - aiptC2 ^ TZTz — = dnHHs 2 - 01/7^4 ■ - — (72) 

02 + a3HH3_4 02 + a3HH3_2 

with some positive constants ai^2,3- Because this is an increasing function of HH.^., and decreasing with 
ptc, it follows that HH3,4 > HH3,2- Rewriting ( [72] ) 

iUis,4{dH - aiptC2 ■ — ^-7777 — ) = HH3,2(ii// - aiptc^ ■ — \— — ) 

a2 + a3HH3 4 02 + a3HH3 o 



and comparing with the Patched values from (66 1, 

HH3^4 dn - aoPTC4,4 ^ ^ 
HH3,2 dn - aoPTC2,2 

for an appropriate positive constant a^. This last inequality shows that PTC2,2 > PTC4^4. Finally, retracing 
back to (64 1, it is not difficult to see that 

ptC2 > ptc^ PTC2,T > PTC4,T- 



Distinct ptC2,ptc^ does not increase robustness On the whole, there are four possibilities to consider: (i) 
PTC2,4 > KpTcci; (ii) PTC2 > KpTcci > PTC4; (iii) PTC4 > Kptcci > PTC2; and (iv) Kptcci > PTC2,4- As 
already mentioned, PTCi > Kptcci in all four situations. 

Situation (i) is similar to the case T2 = T4 already studied, where CI and CN have the form ( 14 1. In case 
(ii), the Cubitus proteins have the form: 

CIl,2 = C/l,2T^, CNi,2 = f/l,2l^ 

CI4 = C/4, CN4 = 0. 

The conditions for wg activation by CI in the second cell require U2 > U4, so the parameters Kch^-g, Hcm-g 
may take values only from components Gc.i or Gc,u, or Gamo- In case (iii) the Cubitus repressor protein is 
not present in the second cell: 

CIi,4 = Ui,4 ^^HaCci ' = ^1.4 l+gagci 

CI2 = U2, CN2 = 0. 

Since CN2 = 0, repression of engrailed on the second cell must now be due to insufficient Wingless activa- 
tion, implying: 

EWG2 < KwG„, < EWG3 

which is impossible, since it was shown that EWG2 > EWG3 for any choice of parameters (see Ap- 
pendixjEj). Finally, in case (iv). Cubitus repressor protein is not present in either the second or fourth cells: 

Cl2,4 = f/2,4, CN2,4 = 0. 

Again it is not difficult to see that Ka,„g, Kcm-g may take values only from components Gc,i or Gen, or Gj^^o 
(it is always necessary that U2 > U4). 

Component Gamo will never become connected to any of the other four, due to opposite requirements on 
the second cell (compare equations (16 1 and ( 17 1). But the comparison above for Clj and CNj show that the 
more general case T2 / T4 only contributes to connect components Gc,i and Gc,u, all the others remaining 
disconnected. 
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G Computing the cylindrical algebraic decomposition 



A CAD for the parameter space G can be computed from equations ([53])-([54)), by imposing the conditions 



h{x) S 3^*^, as given by ( 11 1. By Theorem[T| given h{x) we can solve for EN, EWG, IWG, PTC, and HH 
uniquely as a function of en, wg, ptc, hh and the parameters p. 

First, note that a CAD is not unique, and here we will start by arbitrarily chosing the maximal levels for 
wg, ci, and ptc, that is: 



"owg, aiwG,„g, Ui, U2, U4, T2, Ti (withTi <T2), 



(73) 



with physiological constraints as listed in Table [6] A hierarchy of conditions can then be computed for the 
remaining parameters. 

Second, note that the parameters appearing on the equations for IWG and EWG, as well as those for 
PTC and HH, do not appear on any other equation and, moreover, the unique solution for these four species 
has the same form for any set of parameters (Theorem [TJ. Similarly, all half lives and the Cubitus cleavage 
rate can also be arbitrarily chosen. So, we have a second group of parameters which can be arbitrarily 
chosen, with no conditions to satisfy except for physiological constraints. These parameters are (also listed 
in Table |6]l: 



' exo; 

-f^PTC, -f^HH, [PTC]o, [HH]o, riMPTC, tlmwr, Kvtcw, 



(74) 



Letprf,^^ denote the subfamily of parameters (73 1 and (74i. 

Third, using the computed unique steady state expressions for EN, IWG, EWG PTC, HH, CI, and CN 
write down the conditions for consistency for the expressions of en, wg, ptc, ci, and hh. We have seen that 
EN, CI, and CN have simple expressions: 



and, from Lemma 2.2 



EN^ 

cr^ 

CN*" 



en'^' = (0,0,1,0)' 



(75) 



1 



1 + HciCci 

Hc\Cci 
1 + HciCci 



iUi,U2,0,U^), 

(C/l,C/2,0,C/4). 



The steady state expressions for IWG, EWG, PTC, and HH are more complicated so, for simplicity, we will 
denote them: 

EWG = FEWG = i^EWG(Pnree)> 

IWG = F,wG = i^iwG(Pn,,J, 

PTC = F„c = F„c{Pr,J, 

HH = Fhh = i^HnWeJ. 

We start by showing that there are other parameters which can be arbitrarily chosen, and thus complete the 
proof of Table [61 



Lemma G.l. The parameters Ken^, Kenm,, and Kcnm may take arbitrary values in the interval (0, 1). 
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Proof. The requirements for consistency of the c/*^ expression are: 

EN7 < Ken„ and ENa'' < Ken., and EN*^ > Ken„ and ENj^ < Ken„ 



which are clearly satisfied, in view of (75 1, for any Kenc, £ (0, 1). The requirements for consistency of the 
hh^^ expression are: 



HciCci 



ENr < Ken,,,, [of] Ui—^y^ 
EN2 < Ken,,,, 



> Kc 



or 



1 + HaCc 
ENg"" > Ken,,,, and < Kcn,,,, 



EN4 < Ken,,,, 



or 



Again in view of (75 I, these conditions are all automatically satisfied for any Ken,,,,, Hcmn £ (0, 1) 

Next, the constraints for the parameters in Table [7] are shown. 
Lemma G.2. For system ([T]) with steady state output set ( 1 1 1, the following hold: 

(a) KpTcc, G (0,min{ri,PTC2,T}); 

(b) Kcip,, e (0, ij^HaCa ^i^i^^l' ^^2, Ua}); 

(c) ^cnp.c e ( iSaCc. max{C/i, C/2, C/4}, 1); 

(d) Either KcN» e (0, T^g^ min{[/i, [/a, C/4}) and Kwg., e (0,EWG3), 
or «CN.„ e (0, ifggg^, min{C/i, C/2}) and Kwg» G (EWG4,EWG3). 



Proof. Part (a) follows immediately from Lemma 2.2 since PTC^ = Ti and both VTC12 have to be larger 
than Kpxcci- 

To prove parts (b) and (c), consider the requirements for consistency of the ptc*^'^ expression: 



C/i 



1 



> Kci„,c and C/i 



1 + HaCc 



< K, 



CNpic 



1 + HciCci 

TT ^ A TT ^ClCa 

U2 , , ^ > I^Cl,„c and C/2 < KcNpr. 



1 + HciCr 



' 1 + HriCr 



< Ka,„(. or > KcNptc 

C/4 , , ^ > Kci,«f and Cy4 , , „ ^ < KcNptt- 



1 + HciCc 



1 + HaCc 



The third line is trivially satisfied, while the other lines involve logical ANDs. These immediately yield 
conditions (b) and (c). 

To prove part (d), consider the requirements for consistency of the en^'^ expression: 



EWGi 



EWG2 



< Kwg™ I or I C/l 



< KwG„, I or I c/2 



1 + HciCc 
HciCci 



1 + HciCci 



-C*EWG3 ^ KwGra aud < K^ 



EWG4 



< Kwg™ I or I C/4 



HciCci 
1 + HciCci 
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Theorem|2| says that Fewg4 < ^ewgi = -^ewgs < -^ewg2' so these conditions can be reduced to: 

Kwg™ < i^EWGs and Kcn™ < . , rr mm{Ui , U2 , U4} , (76) 

i + -HciL-ci 



or 



-^EWG4 < ^WGcn < -^EWGa 

and KcNe,, < . , rr mm{Ui,U2}. (77) 

i + -HciL^ci 



It is obvious that the subsets defined by ( 76 1 and ( 77 1 intersect: just choose elements Kwg™ £ (^ewg4 , ^ewgs ) 
and accn™ < i+'£cc, U2, U4}. ' 

Lastly, we come to the parameters in Table [8] and which complete the characterization of the feasible 
parameter space. 

Theorem 4. The set G consists of five disconnected regions of parameters: 

G = Gcj U GcM U GcMI U Gc,/V U Gauio 

each of the regions characterized by Tables^^ |S] 

Proof. The only parameters whose possible intervals have not yet been found are KwGirg> '^cih^;. and Kcm-g- 
Consider now the requirements for consistency of the wg*^ expression. There are three distinct cases, 
depending on wether activation of wg is autocatalytic, or through the CI pathway, or through both. 

Case 1: both CI and IWG contribute to activation of wg. Here wgV = °"..,+«wo.., ^ ^^^^ ^ 



^1 -1 , ^ < I or I Ui > KcN.vg and ^,^^3, < 

I + -HciL'CI ^ + ^ClL-ci / 

TT ^ A TT ^ciCa 

^2 ^ > Kc,,,, and U2 < KcN„, 



< i^chvg I or I > K, 



^CNug ) 



rr ^ ^ I 1 TT ^ciCci 



and 


pCI.WG 
^1WG2 


> '^WGm's 


and 


pCI,WG 
^IWG3 


< '^WGwg 


and 


pCI.WG 


< '^WGwg 



Case 2: only CI contributes to activation of wg Here wg2^ = I'^a'a. ' ^'wg = -^iwc 
1 , I 1 TT HciGi 



Ui , . rr ^ < ^ciMg or ?7i , , ^ > KcN.,.g aud F,wg, < 



1 + HqiCci 1 + -ffciC*, 

1 H G \ 

U2 ^ , „ ^ > Kciwg and ^2 ^ , " " < KcN,.g and F,^^^ < Kwg..^ 

i + -f2ci<-^CI i + -HciOci / 

(0 < KciH.g or I > KcN„g) and Fj^q^ < KwG„.g 



HciGi 



,Tj n < '^"vrg ^4 -I " > KcN„.g ) and F,^^ < KwG„,g- 
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Case 3: only IWG contributes to activation of wg. Here wg2^ = j^^^|r~ , so set Fjwg = ^iwg- 



Ui ^ , ^ < Kci,r« or ^1 , , ^ > KcN.rg and F*q^ < Kwgw 



1 + HciCci 1 + HciCc 



^ ^ ^j^^ HqiCci ^ \ J j-,\vG 

1 + HciCci 1 + HciCci 



U2 , , ^ < Kciu^; or [72 , , ^ > KcNugj and Fj^^^ > ^wcwg 

(0 < Ka„,j or I > KcN,,,,) and F^%.^ < Kwow^ 

rr l , I 1 TT -f^ClC'ci ^ \ 1 r^WG 

^4 ., , „ ^ < Kc,„., , rr ^ > ^cn,,,, and F,^Q^ < KwG,,,. 

i + -nciL'ci J- + -nciL/ci / 



Note that a set of parameters that satisfies case 3 cannot satisfy any of the other two. This is because of the 
conditions on U2- For cases 1 and 2: 

TT ^ ^ J TT HciCci 

while in case 3, the condition to be satisfied is exactly the negation of this. In other words, the region of 
parameter space defined by case 3 cannot be connected to regions defined by cases 1 and 2. 

In addition, note that in cases 1 and 2, U2 = leads to empty intervals for Ka^g, KcNwg)- And a similarly 
conclusion holds when U2 = Ui. This leads to four disconnected regions defined by: Ui > U2 > Ua, 
U2 > Ui, Ui, U2 < Ui, U4 and Ui < U2 < U^. It is now easy to check that, in each of these five regions, 
the intervals for the three parameters Kwcg, Kci„g, and Kcnhs are as given in Table [8] I 

Table 6: Free parameters (physiological constraints, as in f3]). 



Parameter 


Interval 


Ul,2,4, T2 


(0,1] 




(0,^2) 




(0,1) 


Half-lives 


[5, 100] 


{Hx) 




Transfer, cleavage 


(0,1) 


{Cci, Tjxo '"endo, ^M^ f LM ^ 




riMtTc, rMHH, [HH]o, [PTC]o, Kftchh) 




I^ENch ^EN/i/n ^CN/i/i 


(0,1) 
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Figure 3: Projection of set G into the ([/i, U2, U^) space. The regions are defined by the planes U2 = U4 
and U2 = Ui. In region Gauio, Ui can take values in the whole unitary cube, while Ga through Gc.iv do not 
include any of the points in the two planes. In this figure G^mo appears to "intersect" all others, since they 
share values of Ui. However, this is only the projection effect, since not all parameters can be shown. 



Clwg 



CNwg 



Auto 



C, LL 



WGwg 



Figure 4: An example of regions Ge n (solid hne rectangle) and Gamo (dashed line polyhedron). This is the 
projection on the space (Kci„g, {Kcm-g, (Kv/owg), of the fibre over the point represented by "*" in Fig. |3] This 
points corresponds to choosing values for {Ui, U2, U4) in region Gen- 
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kWGen 



kENhh 



H IWG 



aClwg 




kCNhh 



H ptc 



aWGwg 



kWGwg 



kPTCHH 



H PTC 



rendWG 



kClwg 



PTCO 



H ci 



rexoWG 



kCNwg 



HHO 



H CI 



rMWG 



kClptc 



CCI 



H hh 



rLMWG 



kCNptc 



H en 



H HH 



rLMPTC 



kENci 



H EN 



H PH 



rLMHH 



— .'■---^-^M/rt^^-'r>^^VV^-* 



kPTCCI 



H wg 



Figure 5: Parameter histograms out of 70026 parameter sets (refer to model equations for explanation of 
parameters). The notation and scales follow those of Fig. 6 in |[T3l . The half-lifes (denoted Hx) range 
between 5 and 100 mins in a linear scale. The coefficients aCIwg and aWGwg range between 1.0 and 10.0 
also in a linear scale. All other parameters range between 10~^ and 1, in logio scale. 
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QJ 1 0^ ' 0^ ' 0^ ' 0' ' 

50 50 50 50 5 




Parameters 

Figure 6: The distribution functions for the probability of leaving the region (or "mutation"), in each of the 
four regions. is shown in the top row, where the x axis is in logarithmic scale, and t ranges from 

to 40000 (40000 is the maximal number of steps allowed in the random walks). The bottom row shows the 
failed parameters and the percentage of cases where each parameter failed (exit through a certain "face" of 
the polygonal regions). 
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Table 7: Parameters with constraints common to all regions. 



Parameter 



Interval 



'^PTCCI 



(0,min{Ti,Fprc2,T}) 



(0, 



l+HaCa minj^i, ^2,^/4}) 



I^CNptc 



HcNen e (0, T^^;^ min{;7i, U2, Ua}) and Kwg.„ G (0, -Fewgs) 

I^CNeni '^WGen Or 
Kwen g (0, T^^^ min{;7l, U^}) and /twG.n g (-FeWG4) -^EWGs) 



Table 8: The five disconnected components. 



Interval 

(TT HaCa_ I) 



Parameter 



Region 



Gc,I 
U2>U4,Ui 



V^2 i+iJciCci ' ^1 l+iJciCci ^ 



HciCg 



Gc.ii 
Ui>U2> C/4 



(^2i|^,min{f/i,f/4},|gg^^ 



Gc.m 
C/l,C/4>C/2 



I^Wwg 



U4>U2> Ui 



foralH = 1,2,4 



(max{F°,^^3J,l) 



IWGl,2,3,4J 

or 

CI,WG -| pCI,WG\ 
IWG2 ) 



(max{F,wGi,3,4},i^i 



(ma^{FWG^_^j^^wGj 
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